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Abstract 

We develop a numerical scheme to determine which planar snake motions are optimal for 

O" 

CN ' locomotory efficiency, across a wide range of frictional parameter space. For a large coefficient 

U , 

Qh, of transverse friction, we find that retrograde traveling waves are optimal. The optimal snake 
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deflection scales as the -1/4 power of the coefficient of transverse friction, in agreement with an 
asymptotic analysis. At the other extreme, zero coefficient of transverse friction, we propose a 



,^ . triangular direct wave which is optimal. Between these two extremes, a variety of complex, locally 

optimal motions are found. Some of these can be classified as standing waves (or ratcheting 
motions). 
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I. INTRODUCTION 

Snake locomotion has a long history of study by biologists, engineers, and applied math- 
ematicians [l|-l6| . A lack of appendages makes snake motions simpler in some respects than 
those of other locomoting animals |7|. However, snakes can move through a wide range 
of terrestrial p, Isl, |9| , aquatic [10] , and aerial [ll| environments, with different kinematics 
corresponding to the different dynamical laws which apply in each setting. Here we focus on 
terrestrial locomotion, with motions restricted to two dimensions, and a Coulomb frictional 
force acting on the snake. Even with these simplifying assumptions, there are many possible 
snake kinematics to consider. One way to organize our understanding of a diversity of lo- 
comotory kinematics is to propose a measure of locomotory performance (here, efficiency), 
and determine which motions are optimal, and how their performance varies with physical 
parameters. Well-known examples are optimization studies of organisms moving in low- 
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- llTJ I and high- Reynolds-number fluid flows 18|-|23| 



In this work, we use a recently-proposed model for snake locomotion in the plane which 
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has shown good agreement with biological snakes [2|, |5|]. The snake is a slender body whose 
curvature changes as a prescribed function of arc-length distance along the backbone, and 
time. Its net rotation and translation are then determined by Newton's laws, with Coulomb 
friction acting between the snake and the ground. This model is perhaps one of the simplest 
which can represent arbitrary planar snake motions. Previous studies have used the model 
;o find optimal snake motions when the motions are restricted to sinusoidal traveling waves 



5|, or the bodies are composed of two or three links [24]. These motions are represented 
by a small number of parameters (2-10); at the upper end of this range, the optimization 
problem is difficult to solve using current methods. In [1[, Guo and Mahadevan developed 
a model which included the internal elasticity and viscosity of the snake, and studied the 
effects of these and other physical parameters on locomotion for a prescribed sinusoidal 
shape and sinusoidal and square-wave internal bending moments. In the present work, we 
neglect the internal mechanics of the snake, and consider only the work done by the snake 
on its environment, the ground. However, we consider a more general class of snake motions. 
Hu and Shelley assumed a sinusoidal traveling-wave snake shape, and computed the snake 
speed and locomotory efficiency across the two-parameter space of traveling- wave amplitude 
and wavelength for two values of frictional coefficients [5;]. They also considered the ability 



of the snake to redistribute its weight during locomotion. Jing and Alben considered the 
locomotion of two- and three-link snakes, and found analytical and numerical results for the 
scaling of snake speed and efficiency [2^. Among the results were the optimal temporal 
function for actuating the internal angles between the links, expressed as Fourier series with 
one and two frequencies. 

Here we address the optimization problem for planar snakes more generally, by using a 
much larger number of parameters (45 in most cases) to represent the snake's curvature as 
a function of arc length and time. Although our snake shape space is limited, it is large 
enough to represent a wide range of shapes. To solve the optimization problem in a 45- 
dimensional shape space, and simultaneously across a large region of the two-dimensional 
friction parameter space, we develop an efficient numerical method based on a quasi-Newton 
method. The simulations show clearly that two types of traveling-wave motions, retrograde 
and direct waves, are optimal when the coefficient of friction transverse to the snake is large 
and small, respectively. An asymptotic analysis in [25| shows how the snake's traveling 
wave amplitude should scale with the transverse coefficient of friction when it is large. The 
analysis gives a —1/4 power law for the amplitude which agrees well with the numerical 
results shown here. When the transverse coefficient of friction is zero, we give another 
optimally-efficient motion here which is a direct wave of deflection along the body. Between 
these extremes, the numerics show a more complicated set of optima which include standing 
wave (or ratcheting) motions, among a wide variety of other locally-optimal motions. Taken 
together, these results show which planar snake motions are optimal within a large space of 
possible motions and frictional coefficients. 



II. MODEL 
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We use the same frictional snake model as [2|, |5|, |2^ , so we only summarize it here. The 
snake's position is given by X(s,t) = {x{s,t),y{s,t)), a planar curve which is parametrized 
by arc length s and varies with time t. A schematic diagram is shown in figure [1] 

The unit vectors tangent and normal to the curve are s and fi respectively. The tangent 
angle and curvature are denoted 6{s,t) and K,{s,t), and satisfy dgX = cos9, dsV = sin^, and 
K = dgO. We consider the problem of prescribing the curvature of the snake as a function of 
time, K,{s,t), in order to obtain efficient locomotion. When K,{s,t) is prescribed, the tangent 
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FIG. 1: Schematic diagram of snake position, param.etrized by arc length s (nondimensionaUzed by snake 
length), at an instant in time. The tangent angle and unit vectors tangent and normal to the curve at 
a point are labeled. Vectors representing forward, backward and transverse velocities are shown with the 
corresponding friction coefficients /i/, fib, and fit- 



angle and position are obtained by integration: 

9{s,t) = eo{t)+ / i^{s',t)ds', 



x{s,t) = Xo(t) + / cos6{s',t)ds' 



y{s,t) = yo{t) + / sinOis', t)ds'. 
Jo 



(1) 
(2) 
(3) 



The trailing-edge position (xo,yo) and tangent angle 9o are determined by the force and 
torque balance for the snake, i.e. Newton's second law: 



L i-L 

pdttxds = / fr^ds, 
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L pL 

pdttyds = / fyds, 
Jo 
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pX^ ■ dtt^S^ds = / X^-fds. 
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(4) 
(5) 
(6) 



Here p is the snake's mass per unit length and L is the snake length. The snake is locally 
inextensible, and p and L are constant in time, f is the force per unit length on the snake 
due to Coulomb friction with the ground [2|: 

f (s, t) = -pgpt {dtk -njn-pg [pfH{d^ ■ s) + p,{l - H{dtX ■ s))) (dtX ■ s) s. (7) 

Here H is the Heaviside function and the hats denote normalized vectors. When ||i9iX|| = 
we define S^X to be 0. According to ([7j) the snake experiences friction with different 
coefficients for motions in different directions. The frictional coefficients are pf, pb, and pt 



for motions in the forward (s), backward (— s), and transverse (i.e. normal) directions (±n), 
respectively. In general the snake velocity at a given point has both tangential and normal 
components, and the frictional force density has components acting in each direction. A 
similar decomposition of force into directional components occurs for viscous fluid forces on 
slender bodies [26 1. 



We assume that the snake curvature K,{s,t) is a prescribed function of s and t that is 
periodic in t with period T. Many of the motions commonly observed in real snakes are 
essentially periodic in time [2|. We nondimensionalize equations (I4])-(|6]) by dividing lengths 
by the snake length L, time by T, and mass by pL. Dividing both sides by pjg we obtain: 

— / dttxds = / f^ds, (8) 

X^-duXds= / X^-fds. (10) 
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In (l8|)- (fT0!) and from now on, all variables are dimensionless. For most of the snake motions 
observed in nature, L/fifgT"^ ^ 1 [2], which means that the snake's inertia is negligible. 
By setting this parameter to zero we simplify the problem considerably while maintaining a 
good representation of real snakes. fl51)- flTIl become: 

(11) 
(12) 
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f(s,t)= ^' (dtX-h)h (HidtX- 
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(14) 

The equations flTTl) - flT3|) thus involve only two parameters, which are ratios of the friction 
coefficients. From now on, for simplicity, we refer to pt/pf as pt and Pb/f^f as /i^. Without 
loss of generality, we assume /i;, > 1. This amounts to defining the backward direction as 
that with the higher of the tangential frictional coefficients, when they are unequal, pt niay 
assume any nonnegative value. The same model was used in ^ |5|, 12^, and was found to 
agree well with the motions of biological snakes in 



Given the curvature K{s,t), we solve the three nonhnear equations ( !TT|) -( !T3|) at each time 
t for the three unknowns xo(t), yo{t) and 6o{t). Then we obtain the snake's position as a 
function of time by ([I])-(IS])- The distance traveled by the snake's center of mass over one 
period is 

d=\iU <s^ 1) - <s^ 0) ds\ + (j y{s, 1) - y{s, 0) d^ . (15) 

The work done by the snake against friction over one period is 

W= f f i{s,t)-dtX{s,t)dsdt. (16) 

Jo Jo 

We define the cost of locomotion as 

W 

and our objective is to find K,{s,t) which minimize t] at {fib, fit) values which range widely 
over [1, oo) x [0, oo). In our computations, instead of rj, we minimize a different function 
which is more convenient for physical and numerical reasons: 

p = ^g2cos(e(0,l)-e(0,0))_ (^^g^ 

To obtain F from rj we have substituted —d/W for W/d, to avoid infinities in the objective 
function for a common class of motions with d — )■ and W of order 1. The exponential term 
in ( TT8l) penalizes rotations over a period, to ensure that the snake moves in a straight path 
rather than a circular path over many periods. In the Supplementary Material section El we 
give a more detailed explanation for the choice of F. 



III. NUMERICAL MINIMIZATION METHOD 

Our goal is to determine the snake shape, K{s,t), which minimizes F for a given {fiiy,fit)- 
Since K,{s,t) has t-period one, we represent it via a double-series expansion: 

rrai — 1 rii — l 

K{s,t) = ^ '^ {ajk cos{27r jt) + l3jksm{27rjt))Tk{s). (19) 

j=0 k=0 

Here T^ is a Chebyshev polynomial in s: 

Tfc(s) = cos(A;arccos(s)). (20) 

The expansion (fT9|) converges as mi,ni — )■ oo for a large class of n which includes differ- 
entiable functions, and with a convergence rate that increases rapidly with the number of 



bounded derivatives of k [27|. It is reasonable to hope that minimizing F over a class of 
functions flT9|) with small mi and rii gives a good approximation to the minimizers with 
infinite mi and rii. Since the terms in (1191) with coefficients /3ofc are automatically zero, the 
functions given by flT^ are a {2m,i — l)ni-dimensional space. In most of this work we use 
mi = rii = 5, so we are minimizing F in a 45-dimensional space. We give a few results 
at large fit with m,i = ni = 10 for comparison, and find similar minimizers in this region 
of parameter space. The minimization algorithm converges in a smaller portion of {fib, ft-t) 
space as m,i and ni increase, and convergence is considerably slowed by the increased dimen- 
sionality of the search space. We find that irti = ni = 5 is small enough to allow convergence 
to minima in a large portion of (/i;,, fit) space, but large enough to approximate the gross 
features of a wide range of K,{s,t). Another reason to expect that good approximations to 
minimizers can be obtained with small ni is that the dynamical equations flTTl) - flT3l) depend 
only on spatial integrals of body position and velocity. 

Using the dynamical equations (fTTl)-( fT3l) . we show in Supplementary Material section [B] 
that d, W, and F are the same for any reparametrization of time. Such a reparametrization 
can be used to reduce the high-frequency components of a motion while keeping the efficiency 
the same. Therefore, it is reasonable to expect that a good approximation to any minimizer 
of F can be obtained with low temporal frequencies, that is, with an expansion in the form 
of flT9l) with small mi. 



We use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm 23] to minimize F 
numerically, by computing a sequence of iterates which converge to a local minimizer of 
F. The algorithm requires routines for computing the values and gradients of F as well as 
a line search method for finding the next iterate at each step, using the search direction 
provided within the BFGS algorithm. The BFGS algorithm is a quasi-Newton algorithm, 
and therefore forms an approximation to the Hessian matrix of second derivatives, without 
the expense and difficulty of computing the Hessian matrix explicitly. 

We have developed an efficient method for evaluating F and its gradient with respect 
to the {2mi — l)ni shape parameters with just a single simulation of the snake trajec- 
tory over one period. Computing the value of F requires computing the trajectory of the 
snake over one period. We discretize the period interval uniformly with m time points: 
{0, 1/m, 2/m, . . . , 1 — 1/m}. At each time point we write flTT|) - flT3|) as a nonlinear algebraic 
system of three equations in three unknowns: {dtXo,dtyo,dt9o}- To formulate the system. 



we set {xq, yo, Oq} to {0, 0, 0} at each time point, which means that we set the snake in an ar- 
tificial frame rotated from the lab frame. dtOo is the same in the two frames, so after solving 
for {dtXo,dtyo,dt9o} in the artificial frame, we integrate dt9o in time to obtain 6'o(t) in the 
lab frame. This function is the angle between the artificial frame and the lab frame at each 
time, and we use it to rotate {dtXo,dtyo} to the lab frame. We then integrate {dtXo,dtyo} 
in the lab frame to obtain {xo(t),yo{t)}, and integrate K(s,t) using {xo(t) , yo(t) , Ooit)} to 
obtain the body trajectory a posteriori Using {dtXQ,dtyo,dtOo} as the unknowns instead 
of {xo,yo,6o} has two advantages: we avoid cancellation error associated with computing 
discrete time derivatives of {xo,yo,9o}, and we solve m decoupled systems of three equa- 
tions in three unknowns instead of a large coupled system of 3m equations in 3m unknowns 
({xo,i/0)^o} at all times points), which is more expensive to solve. 

We solve the system flTT|) - flT3|) using Newton's method with a finite-difference Jacobian 
matrix. We solve it at each time sequentially from t = 0, using a random initial guess for 
{dtXo, dtyo, dtOo} at t = and an extrapolation from previous time points as an initial guess 
at other t. An important element of our method is the quadrature used to compute the 
integrals in (ITT1)-( IT3|) using ( TT4|) . Here we introduce notation for the tangential and normal 
components of the snake velocity: 

Us = 9tX ■ s ; Un = dt'X ■ h. (21) 

In flT^ we can write 

^■s = ^^^ ; dtX-n=^^^. (22) 

V "2 + ul y/ul + Ul 

These terms have unbounded derivatives with respect to Ug and m„ when Us,Un — )■ 0. The 
local error when integrating ( IT^ using the trapezoidal rule (for example) with uniform 
mesh size h is O (/i^max(us, u„)~2). To obtain convergence in this case it is necessary to 
have h = o(max('Us,'u„)2/^), so the trapezoidal rule (and other classical quadrature rules) 
needs to be locally adaptive when {us,Un) becomes small. 

We use a different approach which allows for a uniform mesh even when {us,Un) — >■ 0. 
We perform the quadrature analytically on subintervals using locally linear approximations 
for Us and u„. Using {dtXo,dtyo,dtOQ} and K{s,t), we compute c^^X, X, and then Us and Un 
on a uniform mesh. On each subinterval, we approximate Us and m„ as linear functions of 
the form As + B using their values at the interval endpoints. We approximate X-*-, s, and 

8 



fi as constants using their values at the midpoints of the intervals. Then the integrals in 
dnD-dnD, using (dl]), are of the form 

Co f ^' + ^ ds (23) 

on each subinterval. We evaluate such integrals analytically using 

h= I . ., ^ = = - log (a2 + 2V^) + log (2 + a2 + 2^1 + 03 + as) , (24) 
Jo ys^ + 02^ + 03 V / 

f sds , /- a^ ^ ,^^s 

h= / , , - = -v/% + VI + Q2 + Q3 - ^h. (25) 

io V s + a2S + 03 2 

When used to evaluate fl23l) . the logarithms in f l24|) -f l25|) have canceling singularities when 
the linear approximations to Ug and Un are proportional (including the case in which one 
is zero), i.e. when AD — BC -^ 0. When this occurs, we use different formulae because 
the integrands are constants or sign functions, which are easy to integrate analytically. 
This quadrature method is 0{h'^) for mesh size h, due to the use of linear and midpoint 
approximations of the functions in the integrands. 

With this method for evaluating (TTT!) - (IT3|) . we compute {dtXQ.dtyo.dtOo] by Newton's 
method as described above, and obtain the snake trajectory and d, W, and F for the 
prescribed curvature function. 

The BFGS method also requires the gradient of F with respect to the {2mi — l)ni 
curvature coefficients {(ypg,Ppq}. For notation, let a be a {2mi — l)ni x 1 vector whose 
entries are the curvature coefficients. We compute a finite-difference gradient using 

_^ F(a+10-®||a||e,) -F(a)) , , 

VFj = ^ J[,\\ / ^^, J = 1, . . . , (2mi - l)ni, (26) 

where e^ is the unit basis vector in the jth coordinate direction. (1261) requires {2mi — l)ni + 1 
evaluations of F, or (2mi — l)ni + 1 computations of the snake motion over a period. We 
employ a very accurate approximation which requires only a single computation of the snake 
motion over a period. 

Computing F(a+ 10^^||a||ej) requires solving b = 0, i.e. ( iTT]) - (iT3l) . for {dtXo,dtyo,dtOo} 
at each time tk = k/m G [0, 1/m, ... 1 — 1/m] using the vector of curvature coefficients 
a+10^^||a||ej. b is a function of {dtXo, dtyo, dtO^Y and the curvature coefficients, so we write 
it ash{{dtXQ,dtyo,dt9QY ,ai+lQ^^\\ai\\ej) and solve b(((9iXo, 9i?/o; <9i6'o)^, a+ 10^®||a||ej) = 
for (c?iXo, dtyo, dtOoY at each time tk- The value of {dtXo, dtyo, dtOoY which satisifies these 



equations at each tk may be written Cfc(a + 10~^||a||ej). Having simulated the snake motion 
and evaluated F when the curvature coefficient vector is a, we have Cfc(a) for each tk- 
Cfc(a) solves h{{dtXo, dti/o, dtOo)^ , a) = 0, a set of equations whose terms are extremely close 
to (within about 10^® of) those of h{{dtXQ,dtyQ.,dtdQ)^ ^s. + 10^^||a||ej) = at each t^. 
Thus c,fc(a) is an extremely good initial guess for Cfc(a + 10~®||a||ej), and we need only 
perform a single Newton iteration with this initial guess to obtain Cfc(a + 10~^||a||ej) to 
machine accuracy (about 10~^^). The one-step Newton iteration is as follows. We define 
Cj = 10^®||a||ej for our description. We start with the usual Taylor series expansion from 
which Newton's method is obtained: 

= b(cfc(a + ej),a + ej) = b(cfc(a), a + Cj) +^(cfe(a), a + Cj) (cfc(a + ej) - Cfc(a)) + 0(||ejf ). 

(27) 
Here J is the 3-by-3 Jacobian matrix of partial derivatives of b with respect to its leading 
3-by-l vector argument: 

Z(x, y),, = 9,^,6,(x, y), z = 1, 2, 3, j = 1, 2, 3. (28) 

Having previously simulated the snake motion with the curvature coefficient vector a to 
evaluate -F(a), we have computed Cfc(a) via a separate instance of Newton's method, and 
we have also calculated and stored J(c,fc(a), a) at each tk- We therefore modify J in (!27|) to 
use this information, while retaining the same order of accuracy in ey. 

= b(cfe(a+ej),a + ej) = b(cfc(a), a + ej)+/(cfc(a), a) (cfc(a + Cj) - Cfe(a)) + 0(||ejf ). (29) 

We solve ( l29l) to obtain Cj^, our approximation to Cfc(a-l-ej) using a single Newton iteration: 

Cjk = Cfe(a) - ^(cfc(a), a)~^b(cfc(a), a + ej) = Cfc(a + ej) + 0(||ej|H. (30) 

( l30l) requires only a single computation of the snake motion over a period, which gives Cfc(a) 
and J(cfc((a), a)^^ at each tk- It also requires (2mi — l)ni evaluations of the function b. 
These are much less expensive than even a single computation of the snake motion over a 
period. The 0(||ej|p) error in Cjk is of the same order of magnitude as machine precision. The 
value of the cost function Fj computed from Cjk is approximately within machine precision 
of F{a + ej) (which corresponds to Cfc(a + ej)), and therefore Fj contributes an error of the 
same order as the ~ 10^^ finite-difference error in (I26p when it is used in place of F(a + ej) 
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in ( 126|) . This level of error leads to essentially no deterioriation in the convergence of BFGS 
up to machine precision, relative to the version with the exact gradient |29]. Therefore we 
obtain the gradient of F for approximately the cost of computing F. Each computation of 
a search direction in BFGS requires only one gradient computation. These ingredients give 
an efficient optimization method, which is critical to finding optima in the high- dimensional 
space of curvature coefficients. 

IV. LARGE-^i RESULTS 
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FIG. 2: Result of one snake optimization run for /if, = 1 and /it — 30. (a) Map of curvature versus arc 
length and time, (b) Values of the objective function (dashed line) and logj^g of the 2-norm of its gradient 
(solid line) with respect to the the iteration count in the optimization routine, (c) Snake trajectory over 
one period. 

In Figure [2] we show the result of one run of the optimization routine with ^^ = 1 and 
lit = 30, starting with a randomly- chosen initial vector of curvature coefficients. In panel 
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a we plot a map of curvature values k(s, t) over one period, and find that the curvature 
approximates a sinusoidal traveling wave with wavelength approximately equal to the snake 
body length. In this case the waves appear to travel with nearly constant speed, shown by 
the straightness of the curvature bands. The bands deviate from straightness near the ends, 
where the largest curvatures occur. In panel b we plot the objective function and log^^o of 
the 2-norm of its gradient over the iteration count of the optimization routine. We find that 
the gradient norm decreases by four orders of magnitude within the first 100 iterations, and 
the objective function reaches a plateau at the same time. In panel c we give snapshots of 
the position of the snake X(s, t) at equal instants in time, moving from left to right, and see 
that the snake approximately follows a sinusoidal path. 
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FIG. 3: Results of three additional snake optimization routines for ^j, = 1 and /ij = 30. (a-c) Maps of 
curvature versus arc length and time after convergence, starting from three different randomly-chosen initial 
iterates. 

Figure [3] shows three different converged results for the same parameters fi^ = 1 and 
fit = 30 but with different randomly-chosen initial iterates. Again we have traveling waves 
of curvature shown by bands, but unlike in Figure [2^, the bands are not straight. Also, the 
bands travel from right to left instead of left to right. The bands are not straight because 
the traveling wave speed changes over time. However, at any given time, the distances in 
s from the maximum to the minimum of curvature are nearly the same in Figure (S^-c and 
Figure [2^, and so are the shapes. In all four cases, the snake passes through the same 
sequence of shapes, but at different speeds, and different (integer) numbers of times in one 
period. Thus the four motions are essentially the same after a reparametrization of time. 
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Correspondingly, the values of F for all four are in the range —5.4418 ±0.0024. The distances 
traveled, d, are within 0.2% of integer multiples of 0.597, with the multiples 4,3,2, and 1 for 
the motions in Figure |2t and Figures |3^-c, respectively, which correspond to the number of 
times they move through the sequence of shapes represented by Figure [3t. The values of F 
are the same whether the waves move from left to right or right to left in s, because /i^ (or 
fih/f^f) = 1, so there is a symmetry between forward and backward motions. 

Over an ensemble of 20 random initializations, this sinusoidal motion yields the lowest 
value of F among the local optima found. We have repeated the search with twice as many 
spatial and temporal modes (mi = 10, rii = 10) and 48 random initializations. About 
three-quarters of the searches converge, and the majority of the converged states are very 
close to those shown in Figures |5] and El The values of F for these states lie in the range 
—5.577 ± 0.026. The mean is systematically lower than the mean for mi = 5, ni = 5 by 
about 2.5%. The decrease in the objective function is slight, and the optimal shape is fairly 
smooth, which indicates that the mi = 5, ni = 5 minimizer is representative of the minimizer 
for much larger mi and rii, at least for /if, = 1 and fit = 30. With mi = 10 and ni = 10, the 
curvature bands of Figures ^ and |3] become straighter near the ends, showing the effect of 
the truncation in the Chebyshev series. However, it is more difficult to obtain convergence 
with four times as many modes, because of the higher dimensionality of the space. Since 
a primary goal here is to obtain the optima across as broad a range of {fit, fJ't) space as 
possible, we proceed with mi = 5 and rii = 5. 

Performing searches with /Xfe = 1 and a wide range of fit, we find that similar sinusoidal 
traveling- wave optima occur for fit ^ 6. In Figure H^ we plot the values of F corresponding 
to the optima over a range of fit on a logarithmic scale, and find that the optima are more 
efficient at higher fit- In panel b, the curvatures of the optimal shapes are plotted for 
the same fit, at the times when the curvatures have a local maximum at s = 1/2. The 
corresponding snake body positions are shown in panel c, and are vertically displaced to 
make different bodies easier to distinguish. In all cases the bodies approximate sinusoidal 
shapes with slightly more than one wavelength of curvature on the body. The curvature 
amplitudes decrease monotonically with increasing fit. 

Performing searches with a wide range of fit and fi^ now, we find that similar sinusoidal 
traveling- wave optima are found for fit ^ 6 and all fib used {1 < fib < 3000). For these 
motions, the curvature traveling waves move backward in s, and the tangential motion is 

13 



a) -3 

-4 

-6 

-7 




C) 



567 10 20 30 60 100 200 300 



b) lOr 



h 





FIG. 4: Traveling- wave optima for ten different fit on a logarithmic scale, a) Values of F. b) Curvature 
versus arc length at the instant when a curvature maximum crosses the midpoint, for the ten fit values in 
(a), c) Snake shapes corresponding to (b). 

purely forward, so they have the same F at all /i;,. However, for /i^ > 1, a separate class 
of (local, not global) optima are also found. These have curvature traveling waves which 
move forward in s, and for which the tangential motion is purely backward. Since friction is 
higher in the backward direction, these local optima have higher F (are less efficient) than 
the oppositely-moving optima, and the values of F for these optima increase with increasing 
/ife, for the same reason. 

In figure [5^ we plot the values of F for these "Reverse" sinusoidal motions at various fij, 
and fit, with leftward-pointing triangles. The triangles connected by a given solid line are 
data for various fit at a particular value of fib, labeled to the right of the line. When fib = 1, 
the reverse sinusoidal motions have the same F values as the forward sinusoidal motions 
shown in figure HI These have the lowest F values across all fib for fit ^ 6, and are thus the 
most efficient motions we have found. 

Three other sets of symbols are used to denote other classes of local optima found by our 
method. All are less efficient than the forward sinusoidal motions, but we present examples 
of them to show some of the types of local optima that exist for this problem. For these three 
sets, shown by circles, downward pointing triangles, and plus signs, we do not join the optima 
with a given fib by lines, to reduce visual clutter. The circles, labeled "Curl," give F values for 
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FIG. 5: Four categories of locally optimal motions: reverse sinusoidal motion, sinusoidal motion with curls, 
reverse sinusoidal motion with curls, and other motions, for ten different [it on a logarithmic scale, a) Values 
of F . Corresponding values of /if, are labeled only for the reverse sinusoidal motions, b) Snapshots of a 
reverse sinusoidal motion with curls. Here /if, = 1 and [ii — 10. c,d) Snapshots of motions in the "Other" 
category, with /i;, = 10 and /it = 60 for (c) and /if, = 6 and /it = 20 for (d). 

shapes which are hke sinusoidal travehng waves, but they have a curl — the snake body winds 
around itself once in a closed loop. Snapshots of the snake motion for such a case are shown 
in figure |5]d. Here the snake moves horizontally, but the snapshots are displaced vertically to 
make them distinguishable. The arrow shows the direction of motion. The tail of the snake is 
marked by a circle. The snake body moves rightward while the curvature wave and curl both 
move leftward. Because the snake overlaps itself, this motion is not physically valid in 2D. 
However, a large variety of such optima are found numerically, with various small-integer 
numbers of curls (and with both directions of rotation). The efficiency of such motions 
decreases as the number of curls increases. The downward-pointing triangles, labeled "Curl, 
Reverse," are for traveling waves with curls which move in the forward direction along the 
snake. These have a yet smaller efficiency for the same reason that the reverse sinusoidal 
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motions are less efficient than the forward sinusoidal motions, in the absence of curls. The 
plusses, labeled "Other," are for other local optima, which have a variety of motions. Panel 
c shows snapshots of a curling (but not self-intersecting) motion in which the snake moves 
horizontally leftward. The horizontal displacement is increased by 8/9 of a snake length 
between snapshots to make them distinguishable. Panel d shows a different curling motion 
in which the snake moves horizontally rightward, and the snapshots are displaced vertically 
to make them distinguishable. These alternative minima increase greatly in number as /ij 
decreases below 6, and some of them outperform the forward sinusoidal motions at these 
transitional fit values. Before discussing the regime of fit below 6, we describe some scaling 
properties for the traveling-wave solution to the problem in the large-fit regime. In a separate 



paper 25| , we derive analytically the form of the optimal prescribed curvature in the limit 
of large fit- In that paper, we calculate that the optimal curvature function is a periodic 
traveling wave with any functional form, but in the limit that the wavelength of the traveling 
wave tends to 0, and such that the snake deflection has a root-mean-square amplitude of 
2^'^ fit ■ The numerically-determined optima in figure H] do not have small wavelengths, 
however, due to the finite number of modes used. This is explained further in Supplementary 
Material section [Cl 

We compare the numerically-found optima to the theoretical results of [251 in figure |6l 
Here we take the plots of figure H] and transform them according to the theory. Figure |6^ 
transforms F from figure UK into 7] by the formula rj = —e'^F~^, since rotations over one 
period are negligible for the traveling- wave solutions. In figure [6|i we plot 77 — 1 versus fit for 
the numerical solutions (squares) from figure UK, along with the small- wavelength theoretical 
result (solid line), which is 

r^^l + ^^;y^ + Oifi-'). (31) 

The numerics follow the same scaling as the theory, but with a consistent upward shift. 
Some degree of upward shift is expected from the finite-mode truncation, which makes the 
numerical result underperform the analytical result. The uniformity of the shift may be due 
to the shape consistently assumed by the numerical optima for various fit. For fit = 300 we 
have also computed rj via the dynamical simulation for a much shorter-wavelength shape 
given by K{s,t) = Ai sm{27f s/Uw + 27rt) with wavelength Uw = 1/8 and Ai chosen so that 
the snake deflection has a root-mean-square amplitude of 2^/^fi~[ . The result is given by 
the open circle, which outperforms the optimum identified in the truncated-mode space, 
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FIG. 6: The /it-scalings of the numerical traveling-wave optima from figure |31 (a) The cost of locomotion 
r] minus unity, on a log scale, for the numerical optima (squares) of figure [3^ together with rj — 1 for the 
analytical optimum (solid line), (b) The curvature from figure |4}d rescaled by /i^ . (c) The deflection from 
figure [It rescaled by /Li/ . 



as expected. The circle is 0.013 above the sohd hne. This deviation is due in part to the 
neglect of the next term in the asymptotic expansion, proportional to yU^^ or 1/300 here. 
Numerically, it is difficult to simulate motions with Uw much smaller than 1/8. Figure |6}3 
shows the curvature, rescaled by /i/ to give a collapse according to [25]. The data collapse 
well compared to the unsealed data in |3)d. Figure Et shows the body shapes from ^ with 
the vertical coordinate rescaled, and the shapes plotted with all centers of mass located at 
the origin. We again find a good collapse, consistent with the curvature collapse. 
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V. MODERATE TO SMALL fif. FROM RETROGRADE WAVES TO STANDING 
WAVES TO DIRECT WAVES 



Near fit = Q there is a qualitative change in the numerically-computed optima. At values 
of /it ^ 6, the global minimizers of F all have backward- moving traveling waves which propel 
the snake forward. In addition, there are local minimizers in which the waves move forward 
along the snake (and the snake moves backward), and overlapping solutions with integer 
numbers of curls appearing in the traveling waves. The values of F at the local minima, 
shown in figure are generally well-separated. In addition, there are a small number of 
local minima with ratcheting types of motions. For fit < Q in figure there are an increasing 
number of local minimizers, denoted "Other," quite close to the global minimizer. These 
motions are perturbed traveling waves. The number of local minima increases greatly as fit 
decreases below 6, and their F values are no longer well-separated. 
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FIG. 7: Values of F for all computed local optima with various random initial iterates, across /ifc at four 
/it (10 (a), 6 (b), 3 (c) and 2 (d)) corresponding to the emergence of perturbed-traveling-wave optima. 

Figure [7] shows the values of F for all computed local optima at four values of fit near 
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the transition: 2, 3, 6, and 10. For each /ij, the F values are plotted against Hb- Also, 
symbols are used to distinguish whether the local maximum in curvature travels forward 
or backward on average, and hence gives the direction of the curvature traveling wave (or 
perturbed traveling wave). At ^t = 10 we find that the smallest F is nearly constant with 
respect to fit,. The same is true for /i^ = 6 except at the largest /if, = 2000, for which there 
are several closely spaced values of F which are slightly below the minimal F at smaller fib- 
These values mark the first appearance of the perturbed traveling waves as fit decreases. For 
fit = 3, the distribution of values is quite different, with a much stronger decrease in F as fib 
increases. There are many clusters of closely-spaced but distinct values of F, showing a large 
increase in the number of local equilibria. At each fib the smallest values of F are spread 
over a band of about 0.5 in width. At fit = 2, the same qualitative trend holds, except that 
the decrease in F with fib is much greater. There is more scatter of values overall and of 
those representing forward-moving traveling waves. There is randomness in the sample of 
F values, particularly evident for fit < 3, because they are values for a finite sample of local 
minimizers starting from random initial k coefficients. 

Intuitively, the minima for fit = 3 represent perturbations of traveling wave motions 
towards "ratcheting" -type motions which involve small amounts of backward motion in a 
portion of the snake to "push" the rest of the snake forward. The appearance of backward 
motion explains the variation of F with fib- As fit decreases from 2000 to 6, F for the 
traveling waves increase, and for fit = 3, other motions outperform the traveling wave 
motions at sufficiently large fib- For fit = 2, the motions tend more toward ratcheting-type 
motions, and we will investigate this further subsequently. 

Another way of viewing the transition is to consider the large- fit solution as fit decreases. 
In the snake motion of figure [St, the snake "slips" tranversely to itself as it slithers forward 
(rightward). This transverse motion has a component in the backward (leftward) direction, 
and the resulting friction with the ground propels the snake forward as a thrust force from 
transverse friction, which balances the drag force from the forward component of tangential 
friction. The thrust force is proportional to fit, and increases as the slope of the snake 
increases, which increases the component of transverse friction in the horizontal direction. 
If fit is decreased, the thrust force can be maintained by increasing the slope of the snake. 
This is indeed observed in the sinusoidal optima as fit decreases, in figure Ht. When fit 
decreases to 3, the slope of the snake becomes nearly vertical at its maximum, and no 
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further decreases in fit can be accommodated by increasing the slope of the snake. Thus, 
near this value of fit, motions different from the sinusoidal traveling wave become preferred. 



^^h 



a) 3000 
1000 

300 

100 

30 

10 

3 

1 



• • • 



••( 






>• • 



0.1 1 10 100 



1 b) 3000 
1000 

300 



0.8 



0.6 




^h 



100 

30 

10 

3 

1 



•• 






' r- m 



0.1 1 10 100 



-1 

-5 
-10 

i" 

-20 
'-25 



FIG. 8: (a) At each (/if,/ib), the maximum value of the function ip over all optima is plotted as a shaded 
dot. (b) The values of F corresponding to the optima in (a) are plotted. 



In figure [8^ we plot a function of the shape dynamics which quantifies the transition away 
from traveling wave motions near fit = 6. The quantity is 

mini maxo.2<s<o.8 I /«(■§, t)| 



^M 



(32) 



maxt maxo.2<s<o.8 \f^is,t)\ 
For K{s,t) which is a traveling wave with a wavelength sufficiently small that a maximum 
or minimum in curvature occurs in 0.2 < s < 0.8 at all times, and for which the maximum 
and minimum are equal in magnitude, ip equals unity. For fit ^ 6, the optimal k{s, t) 
are approximately such traveling waves; we restrict s away from the ends to avoid slight 
intensifications in curvatures near the ends. In figure [8^ we plot the maximum of ip over all 
equilibria located at a given {fit, fib) pair. 

We find ip > 0.8 for fit ^ 6. ip increases slightly from about 0.8 to 0.9 as /it increases from 
6 to 200. For fit = 6, ip decreases slightly but noticeably, from 0.8 to 0.7, as fib exceeds 200, 
showing that at transitional fit, the perturbations away from the traveling wave solutions 
appear first at the largest fib- At fit = 3, the values of tp have decreased to around 0.6 for 
most fib] 'ip is as large as 0.9 at smaller fib and as small as 0.4 at larger fib- For fit < 3, the 
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values of ip are generally much smaller. There is considerable randomness in the values of 
ip plotted for fit < ^ because they are the maxima over samples of local equilibria starting 
from 10-20 random initial k coefficient vectors. 

In figure [H)3 we plot the values of F for the optima whose ^z^- values are shown in panel a. 
We find that F varies smoothly over the region of transitional fit, '^ "£ f^t ^ 10, attaining its 
highest values there. The retrograde traveling-wave motions and other motions are similarly 
efficient at a transitional fit, and the other motions become preferable as fit decreases. 
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FIG. 9: Phase diagram of optimal motions in {fit, fib)-space. 

A global picture of the optimal states is presented in the "phase diagram" of figure [91 
At large fit, the retrograde traveling waves are optimal. At small fit, all of the optima 
have the form of direct traveling waves. Unlike for the large-fit waves, these direct waves 
all have large amplitudes, and their specific form varies with fi^. Near fit = 1 there is 
a region where standing waves with non-zero mean defiection (or ratcheting motions) are 
seen. Examples will be given shortly. These optima co-exist with retrograde-wave optima 
at larger fit. We note that there are other optima which do not fit easily into one of these 
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categories, throughout the phase plane. These other optima are generally found less often 
by the computational search, and are generally suboptimal to the traveling waves, but there 
are exceptions. For the sake of brevity we do not discuss such motions here. 

To obtain a clear picture of the snake kinematics across the transitional region shown in 
figure ini we give figures showing snapshots of snake kinematics, together with spatiotemporal 
curvature plots, at /i;, = 1 and 3 in this section and at fib = 2, 5, and 20 in Supplementary 
Material section [Dl Each of these figures represents one horizontal "slice" through the phase 
diagram in figure [H 
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FIG. 10: Snapshots of optimal snake motions at various fj,t (one per column, labeled at the top), for /i^ = 1. 



In figure [THl we show the optimal snake kinematics for /i^ = 1 and eight fit spanning the 
transition. Each set of snapshots is arranged in a single column, moving downward as time 
moves forward over one period, from top to bottom (shown by arrows). Similarly to figure 
lib and d, the snapshots are shown in a frame in which the net motion is purely horizontal, 
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from left to right. However, the snapshots are given an extra uniform vertical spacing to 
make them visable (i.e. to prevent overlapping). Below the snapshots, the corresponding 
curvature plot is shown, in the same format as in figures [2] and [31 with s increasing from to 
1 on the horizontal axis, from left to right, and t increasing from to 1, from bottom to top. 
We have omitted the axes, and the specific curvature values for each plot, instead using a 
single color bar (shown at bottom) with curvature values which range from the maximum to 
the minimum for each plot. By omitting the axes we are able to show several of the optimal 
motions side-by-side. By comparing a large number of motions in a single figure, it is easier 
to see the general trends. 

We begin by discussing the rightmost column in figure dUl with /ij = 6. A retrograde 
traveling wave is clearly seen, so this optimum is a continuation of the large- /xj solutions. At 
this fit, just above the transition, the wave amplitude is large, and the wave clearly moves 
backwards in the lab frame. Thus the snake slips backward considerably even though it 
moves forward, because fit is not very large. Moving one column to the left, the optimum 
at /it = 3 shows a clear difference. At certain instants, the curvature becomes locally 
intensified near the snake midpoint. This occurs when the local curvature maximum reaches 
the snake midpoint. The motion is still that of a retrograde traveling wave, but now the 
wave shape changes over the period. Hence the displacement is not a traveling wave of 
the form g{s + Uu,t) for a fixed periodic function g. Comparing the curvature plots for 
fit = Q and 3, we see that the bands have nearly constant width across time for fit = 6. 
At /it = 3, by contrast, the curvature bands become modulated. The bands show bulges 
near the curvature intensifications. Moving leftward again, to fit = 2 and then to 1.5, the 
curvature intensifications increase further, showing a general trend. In all cases, however, 
the curvature peaks move backward along the snake, even as the peak values change, so 
these motions are retrograde waves. At fit = 1.2, a novel feature is seen: the endpoints of 
the snake curl up (but do not self- intersect), forming "anchor points." Between these points, 
the remainder of the snake moves forward. In the curvature plot, bands can still be seen 
corresponding to retrograde waves. The three remaining columns, moving leftward, all show 



direct waves. That is, the wave moves in the same direction as the snake 30|]. No optimum 
is shown for fit = 1. The point fib = l,fit = 1, in which friction is equal in all directions, is 
in some sense a singular point for our computations. Near this point, it is difficult to obtain 
convergence to optimal motions, and optima are difficult to classify clearly as retrograde 
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or direct waves (or another simple motion). Nonetheless, we have computed one optimum 
which can be classifed as a direct wave, but it has complicated and special features which 
are not easy to classify. In the three left columns, the direct waves can be recognized by the 
bands in the curvature plots. The bands have the direction of direct waves. However, the 
bands are not nearly as simple in form as those for the retrograde waves, at large /ij. 
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FIG. 11: Snapshots of optimal snake motions at various /it (one per column, labeled at the top), for /i;, = 3. 



Figure [TT] shows optima for /i(, now increased to 3. The same transition is seen, but now 
in two middle columns, with /ij = 1 and 1.5, standing waves are seen. The curvature plots 
are nearly left-right symmetric. In the corresponding motions, the snake bends and unbends 
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itself. Because //;, > 1, there is a preference for motion in the forward direction even though 
the curvature is nearly symmetric with respect to s = 1/2. These optima did not occur 
for yUf, = 1 (figure [TU1) because symmetric motions with equal forward and backward friction 
would yield zero net distance traveled. It is somewhat surprising that such motions can be 
completely ineffective for /i^ = 1 and yet represent optima for /if, as small as 2, as shown in 
figure m 

We have presented only a limited selection of optima for moderate fit, and largely with 
pictures, because it is difficult to describe and classify many of the motions we have observed. 
Future work may consider this interesting region of parameter space further. 

We now move on to the limit of zero fif In this limit, we find a direct- wave motion with 
zero cost of locomotion. This explains the preference for direct waves at smaller fit in this 
section. 

VI. ZERO nt: DIRECT- WAVE LOCOMOTION AT ZERO COST 

When fit = 0, it is possible to define a simple shape dynamics with zero cost of locomotion. 
The shape moves with a nonzero speed while doing zero work. It is a traveling wave with 
the shape of a triangle wave which has zero mean vertical defiection: 

y{s, t) = Ao-l sgn(sin(27r(s - t)))ds. (33) 

The vertical speed corresponding to f l5^ is a square wave: 

dtvis, t) = -AoSgn(sin(27r(s - t))). (34) 

The wave of vertical defiection moves forward with constant speed one, and the horizontal 
motion is also forward, with a constant speed U: 

dtx{s, t) = U (35) 

We determine U by setting the net horizontal force on the snake to zero: 

I dix -83^(13 = 0. (36) 

with 

^"^ I (37) 

74osgn(sin(27r(s — t))) 
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Solving ( 136|) for U we obtain 

U = , ° ■ (38) 

With this choice of U, S^X ■ s is identically zero on the snake, so the motion is purely in the 
normal direction. Therefore, there are no forces and torques on the snake, and no work done 
against friction. However, the snake moves forward at a nonzero speed U. This solution has 
infinite curvature at the peaks and troughs of the triangle wave, and we may ask if zero cost 
of locomotion [i]) can be obtained in the limit of a sequence of smooth curvature functions 
which tend to this singular curvature function. To answer this question, we have set y to 
the Fourier series approximation to (1331) with n wave numbers, for n ranging from 2 to 1000. 
At each time t, we have then calculated a net horizontal speed U{t), net vertical speed V{t) 
and a net rotational speed dR{t)/dt such that Fxif) = 0, Fy{t) = and r{t) = 0. We find 
very clear asymptotic behaviors for U{t), V(t), dR(t)/dt, and rj: 

Uit) = -^=^ + 0{n-') , V{t) = 0{n-') , ^ = 0{n-') , v = 0{n-'). (39) 

Therefore zero cost of locomotion does in fact arise in the limit of a sequence of smooth 
shape dynamics. This type of traveling wave shape dynamics, which yields locomotion in 
the direction of the wave propogation, has been called a "direct" wave in the context of 



crawling by snails [30[ and worms 31|. 



The triangular traveling wave motion is particularly interesting because it can represent 
an optimal motion at both zero and infinite fit, and its motion can be computed analytically 
for all fif We show how this motion embodies many aspects of the general problem in 
Supplementary Material section [El 

VII. CONCLUSION 

We have studied the optimization of planar snake motions for efficiency, using a fairly 
simple model for the motion of snakes using friction. Our numerical optimization is per- 
formed in a space of limited dimension (45), but which is large enough to allow for a wide 
range of motions. The optimal motions have a clear pattern when the coefficient of trans- 
verse friction is large. Our numerical results indicate that a traveling-wave motion of small 
amplitude is optimal, which agrees with an asymptotic analysis in [25;] . When this coefficient 
is small, another traveling-wave motion, of large amplitude, is optimal. When the coefficient 

26 



of transverse friction is neither large nor small, the optimal motions are far from simple, and 
we have only begun to address this regime. 

A natural extension of this work is to include three-dimensional motions of snakes. In 
other words, one would give the snake the ability to lift itself off of the plane [l|, S, |5[, as 
occurs in many familiar motions such as side- winding [9]. Another interesting direction is 



to consider the motions of snakes in the presence of confining walls or barriers 



m. 
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Supplementary material for: 

Optimizing snake locomotion in the plane. I. Computations 

Appendix A: Supplementary Material: Objective function 

To understand our choice of F, we first describe the class of motions X(s,t) which result 
from periodic-in-time K{s,t). Let X(s,t) and (9(X(s,t) be given for the moment. Let us 
rotate X(s,t) by an angle a and translate X(s,t) by a vector v, uniformly over s, so we 
have a rigid body motion. Then s and n are also uniformly rotated by a. Let us also rotate 
(9iX(s,t) uniformly by a. Then by flT4|) . f is uniformly rotated by a. Hence, if X(s,t) and 
(9tX(s,t) are such that the force and torque balances (TTTl) - (TT3l) hold, they will continue to 
hold if we apply a rigid-body rotation and translation to X(s,t) and rotate 9iX(s,t) by 
the same angle. Now, since K,{s,t + 1) = K,{s,t), X(s,t + 1) is obtained from X(s,t) by a 
rigid-body rotation by an angle ai and translation by a vector vi. In fact, the trajectory 
of the body during each period is the same as during the previous period but with a rigid 
translation by vi and rotation by ai. Over many periods such a body follows a circular 
path when ai 7^ 0, and because of the rotation, the distance traveled over n periods may be 
much smaller than n times the distance traveled over one period. To obtain a motion which 
yields a large distance over many periods, we restrict to body motions for which X(s, t + 1) 
is obtained from X(s,t) by a translation only, with rotation angle ai = 0. We implement 
this constraint approximately in our optimization problem, using a penalization term. We 
multiply d in flT7|) by a factor which penalizes rotations over a single period: 

rf^rfe2'=-W0'i)-''(0'0)) (Al) 

flAll) involves the change in angle at the snake's trailing edge over a period, which is the 
same as ai. The exponential factor has a maximum of e^ when ai = 0, and a minimum of 
e~^ when ai = vr. Using the factor of 2 in the exponent of flAll) . all of the optima found 
by our computations involve very little rotation — less than 0.5 degrees over all of (/i6,/if) 
space. When the factor is decreased from 2, some of the optima begin to involve nontrivial 
amounts of rotation. In short, the factor of 2 represents a choice of how much to weight 
lack of rotation relative to distance and work in the objective function. As described in the 
results section, for fit ^ Q and all fib, we find optima which have an analytical traveling wave 



form, and in the limit of large fit, rotations tend to zero for these solutions. These optima 
are found even with no rotation penalty, so the rotation penalty plays a negligible role in this 
region of parameter space. For smaller fit, however, in the absence of the rotation penalty, 
many of the optima have significant rotations. 

A second modification to rj is made to avoid large values of 77. There are snake motions 
(i.e. with the K{s,t) = k(1 — s,t) symmetry and fib = 1) which have d = and W ^ 0, 
giving 1] = 00. The minimization search commonly encounters motions which are nearly 
as ineffective throughout shape space, and t] has a large gradient in the vicinity of such 
motions. It is challenging to minimize such a function numerically. We therefore substitute 
—d/W for W/d. The minima are the same for both because each is an increasing function of 
the other, but —d/W varies much more smoothly in most of shape parameter space. There 
are cases where —d/W tends to —00, but these only occur in a small region of parameter 
space which is rarely approached by our computations. 

Combining these ideas, the function we minimize in place of 77 is 

p = ^g2cos(0(o,i)-e(o,o))_ ^j^2) 

Appendix B: Supplementary Material: Invariance under reparametrization of time 



Here we show that d, W, and F are the same for any reparametrization of time — in other 
words, for k,{s, 7(t)) with any nondecreasing differentiable mapping 7 from [0, 1] to [0, 1] with 
7(0) = and 7(1) = 1. Such a reparametrization can be used to reduce the high-frequency 
components of a motion while keeping the efficiency the same. Therefore, it is reasonable to 
expect that a good approximation to any minimizer of F can be obtained with low temporal 
frequencies, that is, with an expansion in the form of ( TT9l) with small mi. To show that d, 
W, and F are the same for a reparametrization of time, let k{s, t) be given and 7(t) be the 
aforementioned parametrization function. Let X(s, t) be the motion obtained by integrating 
K{s,t) from s = as in (II])-(l3]) with the integration constants Xo(t) = X(0,t) and 6'o(t) 
chosen to satisfy the dynamical equations flTT|) - flT3|) . Now we show that X(s, t) = X(s, 7(t)) 
also satisfies ( fTT]) - (fT3l) at each time t. By the chain rule of differentiation, 

r5 dX{s,',{t))/dt _ 7'(i)a,X(s,7)|,.,(i) _g^ ,| , 



using 7'(t) > 0. Let i{s,t) be the force distribution corresponding to X(s,7(t)). It is given 
by (fT4l) with S^X in place of 9(X and (s, n) in place of (s, n), where (s, n) are the unit vectors 
tangent and normal to X(s,t). Using fIBll) we obtain f(s,t) = i{s,-y(t)). Since flTT|) - flT3|) 
hold for f (s, t) for all t G [0, 1], they also hold for f (s, 7(t)) for all t G [0, 1] since 7(t) G [0, 1]. 
Thus they hold for f(s,t) in place of f(s,t), so X(s,t) satisfies the dynamical equations 
(ITT])-flT3l). The distance traveled d for this motion is 




d=M / dtx{s,t)dsdt] +1 / dty{s,t)dsdt] . (B2) 

(9^a;(s,7)(is(i7 J +( / / f^7l/(s,7)(is(i7 J (B3) 

(B4) 

using flT5|) . The work done is 

W= f f i{s,t)-dtX{s,t)dsdt= f f f{s,-f)-d^X{s,-f)dsd-f = W (B5) 

Jo Jo Jo Jo 

using ( IT6l) . The net rotation over a period is 

^(0,1) -^(0,0)= / dt9i0,t)dt= [ d^e{0,-f)d^ = 9(0,1) -9(0,0), (B6) 



7*- 
^0 



so F (TTSll is unchanged by the reparametrization of time. 

Appendix C: Supplementary Material: Effect of finite mode truncation on large-^j 
optimum 



The theory of [25|] predicts that the wavelengths of the optimal shapes should tend to 
zero, but the wavelengths for the computed optima shown in figure Hb and c are all in the 
neighborhood of one body length, not particularly small. The reason is that it is not possible 
to represent traveling waves with very small wavelengths with the value of rii (5) used here, 
which includes polynomials in s of up to fourth degree. 

A general periodic curvature function can be represented as a Fourier series in x = 
s + Uwt. For any such function, with amplitude scaled appropriately, |25| predicts that the 
efficiency should tend to the optimum as Uw — ;■ 0. For simplicity we focus on a single mode 
sin(27rx) = sm(2ns/Uy^ + 27rt) and ask: how well can it be approximated by the truncated 
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FIG. 12: The approximation error in representing (a) sine and (b) cosine functions with different wavenum- 
bers using either five or ten Chebyshev modes (corresponding curves are labeled) . Panels c-e show snapshots 
of curvature versus arc length in the optimal numerical traveling wave motions with either five (dashed line) 
or ten (solid line) Fourier and Chebyshev modes. The snapshots are at times when there are local curvature 
maxima at (c) s = 0.25, (d) s = 0.5, and (e) s — 0.75. 



Chebyshev series we have used in the numerics? At any time t, sin(27rs/f/^ + 27rt) can be 
decomposed into a hnear combination of sin(27rs/f/i„) and cos(27rs/f/^), so we compute the 
best approximations to these two functions in L°° norm by Chebyshev series with n\ terms. 
In figure [12] we show the resulting approximation errors (in L^ and L°°) for ni = 5 and 10, 
and for sin(27rs/f/^) (panel a) and cos(27rs/f/^) (panel b). We find the error transitions from 
small to large for the wavelength [/^ in the neighborhood of 1 for rii = 5 and 1/2 for n\ = 10. 
For comparison, in panels c-e we plot snapshots of the curvatures from the optimal snake 



motions identified previously for fit = 30 witli mi = rii = 5 and mi = rii = 10. Tlie panels 
show the snapshots at instants when the curvatures have local maxima at s = 0.25 (c), 0.5 
(d) and 0.75 (e). These shapes are approximately sinusoidal, and those with mi = ni = 10 
have a somewhat higher wavelength than those with ttii = rii = 5. The wavelengths for 
mi = ni = 5 appear to be somewhat larger than 1, while those for rrii = ni = 10 appear to 
be approximately 1.5. Based on panels a and b, we would have expected the m,i = ni = 10 
cases in c-e to have somewhat shorter wavelengths, but in the space spanned by the finite 
number of Chebyshev polynomials, the smallest Uw is not the only consideration: certain 
functional forms for the curvature should be expected to perform better than others, and 
the results here may tend towards such forms. 

Appendix D: Supplementary Material: Optimal snake motions across the fit tran- 
sition 

Figure [13] shows optima for /if, = 2, in the same format as figure [T0| at a different set of 
seven fif The motions occur on a horizontal line in figure [9] which touches the bottom of the 
region in which standing waves occur. Starting in the right column, which has /it = 3, we 
see a retrograde wave, with modulations (or bulges) in the bands of the curvature plot. In 
the next five columns, we see the transition from retrograde to direct waves. However, the 
curvature plots begin to resemble a checkerboard pattern. This resemblance is incomplete, 
and the motions are still retrograde or direct waves, but not standing waves. 

Figure [H] shows optima in for /i^ increased to 5. The retrograde wave is somewhat difficult 
to discern at /it = 3, without consulting the curvature plot. The motions are qualitatively 
similar to those for fij, = 3, shown in figure fTTl 

Figure [TS] shows optima in for /i^ increased to 20. For fit = 0.2 we show two local optima. 
The second column shows a "swing-and-anchor" motion. The snake forms a loop at one 
end, the "anchor," and then the rest of the body swings forward, rotating about the anchor. 
The motion also shows, briefiy, a direct wave. This example shows that optima can combine 
multiple types of motions within a single period of motion. We have seen many other such 
"mixed" optima at moderate fit, but they are not reproduced here for brevity. The other 
fit = 0.2 optimum, shown in the first column, is essentially a direct wave. The optima for 
fit = 0.3, 0.7, and 3 also involve the formation of "anchors" near the endpoints, somewhat 
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FIG. 13: Snapshots of optimal snake motions at various yn (one per column, labeled at the top), for /i;, — 2. 



similar to the concertina motion of biological snakes [4]. The optimum at /ij = 2 is a very 
symmetrical ratcheting motion. 
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FIG. 14: Snapshots of optimal snake motions at various /xt (one per column, labeled at the top), for ^b — 5. 



Appendix E: Supplementary Material: Triangular ^vave motion across /ij 



We again consider the triangular wave motion (I33ll- (l35ll . but now /ij is nonzero, so ([36 
takes the more general form 



Hdt^ ■ ssj; + fitdt'^ ■ hn^ds = 0, 

n = (fifH{dtX ■ s) + fi,{l - H{dtX ■ s)) 
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FIG. 15: Snapshots of optimal snake motions at various /it (one per column, labeled at the top), for /i^ — 20. 



with s as in fl371) and 



n 



-AoSgn(sin(27r(s — t))) 



'1-Al 



We solve (EI]) for f/: 



t/ 



n + if^t- n)Ai ■ 



(E3) 



(E4) 



When /it = 0, (lE4p becomes (138|) . The tangential speed is 



Us 



Usx + dtysy 



<0, 



(E5) 



so "H = ^b- The snake speed (1E4|) . as well as Us and m„, are functions of Aq and ^t/fJ^b 
only. By flT6l) and flT7|) . so is ////ib- In figure [T6l we give a contour plot of logiQ?7/yUb as a 
function of g = arcsin Aq and log^g fi't/fJ'b- Q is the angle that the triangle wave makes with 
the horizontal axis. 




FIG. 16: Contour plot of logj^g 77//ib as a function of the triangular wave angle q = arcsin ^o and the friction 
coefficient ratio logj^g ^itl l^b- 



The forward speed U is zero, and the cost of locomotion infinite, when ^t = IJ-b, repre- 
sented by a horizontal line near the middle of figure [161 This is an example of the "singular 
point" mentioned in the previous section, when /^t = A^b = /^/ = 1- When /i^ > fi^, above 
the horizontal line, ^ < 0, and the snake moves as a retrograde wave, rj/^h decays to 1 as 
fJ't/fJ'b -^ oo. In this regime, for fixed fit, V itself is the product of fib with a function that 



decays more slowly than /if, as /i;, — )■ 0. Thus t] is ininiinized when /Xb — )■ 1, the minimal value 
of fib- We plot the optimal amplitude in the limit of large fit/fJ'b (derived for a wave of any 
shape in (25|]) as a dashed line at the upper left of the figure, with the formula for the optimal 
amplitude in this limit. In the lower portion of figure Uni yU* < fib, U > 0, and the snake 
moves as a direct wave. Here r]/fib decays to as fit/fJ'b -^ 0. In this regime, for fixed /i^, 
1] itself is the product of /ib with a function that decays more rapidly than fx^^ as Hb — ?■ C)0. 
Thus rj is minimized by taking /if, — )■ oo and fit -^ 0. The wave angle which minimizes 7] can 
be found by directly minimizing t] with respect to q. We find that in the limit fit/f^b -^ 0, 
the optimal wave angle tends to 7r/2, i.e. the optimal wave becomes steeper. 
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